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Introduction 


Hyperthermia  has  been  shown  to  be  an  efficacious  adjuvant  therapy  in  the  treatment 
of  recurrent  breast  cancer.  [1,2, 3]  One  of  the  most  flexible  and  proven  methods  to  induce 
hyperthermia  is  ultrasound.  In  order  to  optimally  utilize  such  clinical  devices  it  is  essential 
that  the  operator  have  the  ability  to  accurately  plan  a  treatment.  The  clinical  hyperthermia 
situation  is  diagrammed  in  Figure  1.  A  treatment  plan  consists  of  a  defined  motion  path 
and  applied  transducer  power  that  will  raise  the  temperature  of  the  tumor  to  a  desired  set 
point  for  a  specified  time,  while  minimizing  the  effect  on  surrounding  tissue.  The  treatment 
plan  must  be  adjusted  for  tissue  geometry,  absorption,  and  desired  temperature  history 
specific  to  each  patient.  Planning  a  treatment  begins  with  the  accurate  estimation  of  the 
absorbed  power  field  created  by  the  ultrasound  transducer.  Thus,  the  ability  to  accurately 
estimate  this  field  is  critical  to  producing  a  proper  and  effective  treatment  plan. 

Goals 

It  is  the  goal  of  this  research  to  provide  improved  numerical  modeling  of  ultrasound 
propagation  through  breast  tissue  and  the  post-mastectomy  chest  wall  for  use  in  patient 
treatment  planning  of  hyperthermia  cancer  treatments.  With  this  improved  model  the 
clinician  can  plan  an  ultrasound  hyperthermia  treatment  so  as  to  produce  the  required 
thermal  dose  while  minimizing  patient  pain  due  to  excessive  temperature  or  ultrasound- 
tissue  interactions.  A  clinical  system  is  being  constructed  to  combine  model  results  with 
experiments  in  an  inverse  technique  to  estimate  the  absorbed  power  field  in  the  breast  and 
chest  wall  during  ultrasound  hyperthermia.  This  improved  planning  will  aid  the  clinicians 
in  providing  better  hyperthermia  treatments  which  will  improve  treatment  response  rates. 

The  system  under  development  consists  of  the  following  subsystems: 

1.  An  ultrasonic  B-mode  imager  based  clinical  data  acquisition  system  to  obtain 
patient  anatomy,  measure  attenuation,  and  measure  absorbed  power. 

2.  A  clinical  technique  using  the  ultrasound  treatment  system  to  directly  measure 
absorbed  power  at  specific  points. 

3.  An  improved  model  of  ultrasound  propagation  through  breast  tissue. 

4.  An  inverse  technique  that  integrates  experimental  measurements  and  modeling 
results  to  obtain  the  “best"  prediction  of  the  ultrasound  power  deposition  field. 


The  B-mode  diagnostic  ultrasound  imager  in  Subsystem  1  is  used  to  construct  a 
geometric  model  based  upon  each  patient's  anatomy  and  to  measure  the  ultrasound 
attenuation  coefficient  at  selected  areas  within  the  tissue  region.  In  Subsystem  2  the  clinical 
treatment  system  is  used  in  a  scanning  protocol  to  locate  thermocouples  imbedded  in  the 
breast  tissue.  The  response  of  these  thermocouples  to  the  scanned  ultrasound  gives  a  direct 
measure  of  the  power  deposition  at  that  point.  Subsystem  3  is  an  ultrasound  propagation 
and  power  deposition  model.  For  this  model  we  are  currently  considering  2  options:  a) 
Hybrid  Model  -  Ultrasound  propagation  from  the  transducer  into  the  breast  is  modeled  with 
a  hybrid  of  Green's  function  solutions  used  in  the  homogenous  water  region  and  a  finite 
element  solution  of  the  acoustic  wave  equation  in  the  inhomogeneous  breast  region,  b) 
Parabolic  Model  -  This  is  a  forward  marching  parabolic  model  that  is  quite  fast,  but  does 
not  take  into  account  reflected  energy  at  interfaces.  Subsystem  4  is  an  iterative  inverse 
technique  that  optimizes  the  parameters  of  the  power  deposition  model  by  forcing  the  model 
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1  k. 


results  to  match  the  experimental  measurements.  The  model  is  used  to  predict  the  power 
deposition  at  the  thermocouple  locations,  the  error  between  model  and  experiment  is 
measured,  and  the  error  values  are  used  to  adjust  the  model  parameters  for  the  next 
iteration. 


Patient 


Membrane 


Water  Bath 


Ultrasound 

Propagation 


Ultrasound  Transducer 

Figure  1 :  Diagram  of  the  clinical  situation  for  a  hyperthermia  treatment  of  the  breast  or 
chest  wall  area. 
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Revised  Task  Description  and  Milestones: 

As  in  any  research  project,  during  the  course  of  the  project  the  researchers  often 
discover  that  some  tasks  are  easier  than  expected;  other  tasks  are  difficult  or  impossible  and 
are  replace  with  alternative  solutions.  As  can  be  seen  from  the  project  description  in  the 
introduction,  this  has  certainly  occurred  here.  Difficulties  with  the  initial  modeling 
techniques  have  resulted  in  a  search  for  alternatives.  Investigation  of  the  Geometry 
Acquisition/Registration  problem  has  resulted  in  expanded  goals  in  this  area.  In  an  attempt 
to  clarify  what  has  been  accomplished  and  what  is  left  to  accomplish  in  this  research 
program,  the  original  Task  Description  and  Milestones  chart  has  been  modified  to  reflect 
the  current  project  tasks  and  goals. 
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Outline  of  Accomplishments: 

1st  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Acquisition  of  ultrasonic  B-scan  instrument. 

•  Development  of  interface  hardware  to  isolate  analog  A-scan  information. 

•  Acquisition  of  high  speed  A/D  system  and  development  of  control  software. 

•  Investigation  of  reflection  based  acoustic  attenuation  measurement  techniques. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Initial  discussions  and  verification  of  equipment. 

Subsystem  3:  Model  Development 

•  Development  and  testing  of  a  2D  finite  element  based  model. 

•  Identification  of  alternative  solutions. 

Subsystem  4:  Inverse  Technique 

•  Initial  discussions. 
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2nd  Year  Accomplishments: 


Subsystem  1:  Geometry  Acquisition  and  Attenuation  Measurement 

•  Development  of  Log  Spectral  Difference  software  for  measuring  the  acoustic 
attenuation  from  a  single  A-line. 

•  Development  and  testing  of  a  Force/Balance  test  apparatus  for  direct 
measurement  of  acoustic  attenuation. 

Subsystem  2:  Experimental  SAR  Measurement 

•  No  work  accomplished. 

Subsystem  3:  Model  Development 

•  Extension  of  the  2D  FE  model  to  three  dimensions.  Combination  of  the  FE 
model  with  an  integral  method  model  to  improve  the  speed  of  the  model. 

•  Development  of  a  3D  parabolic  model  as  a  much  faster,  possibly  less  accurate 
alternative. 

•  Development  of  a  3D  Green’s  Function  model  for  use  as  the  “gold  standard” 
with  which  the  other  models  will  be  compared. 

•  Initial  model  comparisons. 

Subsystem  4:  Inverse  Technique 

•  No  progress. 


First  Year  Summary 

The  first  year  efforts  centered  on  initial  development  of  the  ultrasound  propagation 
models  and  acquisition  and  initial  development  of  the  diagnostic  ultrasound  equipment.  A 
brief  discussion  of  these  efforts  is  included  below.  For  more  detailed  information  refer  to 
the  first  year  progress  report. 

Initial  modeling  efforts  focused  on  two  potential  methods:  Green’s  Function 
Method  where  the  transducer  face  is  modeled  as  a  distribution  of  point  sources,  and  the 
Finite  Element  Method  where  FE  techniques  are  used  to  solve  the  acoustic  wave  equation. 
A  two  dimensional  FE  code  was  written,  and  verification  runs  were  under  way  at  the  end 
of  the  first  year.  The  results  of  these  runs  and  further  investigation  of  the  modeling 
possibilities  led  to  the  development  of  a  hybrid  model  that  combines  the  features  of  the  two 
techniques  to  provide  improved  speed  and  accuracy.  Details  of  this  model  are  discussed  in 
the  second  year  accomplishments  below. 

During  the  first  year  a  diagnostic  ultrasound  B-mode  imager  was  acquired.  This 
imager  was  modified  to  allow  access  to  the  raw  RF  waveform  (A-line)  for  each  scan  line  of 
the  B-mode  image.  A  schematic  of  the  resulting  system  is  shown  in  Figure  2.  The 
ultrasound  unit  (Scanner  250,  Pie  Medical  USA  3535  Route  66,  Neptune  NJ  07753)  was 
modified  to  output  the  RF  signals  along  with  sync  signals.  A  custom  hardware  selection 
device  was  constructed  to  read  these  signals  and  output  a  single  A-line.  A  PC  based  data 
acquisition  system  was  acquired  that  allows  digitization  of  the  A-line  signals.  Once  the  data 
is  in  the  PC  it  is  available  for  use  in  software  techniques  to  estimate  the  ultrasonic 


attenuation  field  in  the  tissue.  Along  with  the  hardware  development,  software  was 
developed  to  operate  the  hardware  and  perform  the  attenuation  estimate.  [4] 


Figure  2:  Diagram  of  the  B-mode  ultrasound  imager  based  clinical  data  acquisition  system 
used  to  obtain  patient  geometry  and  measure  attenuation. 


Second  Year  Accomplishments 

The  second  year  accomplishments  will  be  discussed  in  terms  of  the  four 
subsystems  identified  above. 

Subsystem  1:  Ultrasound  Based  Anatomy  and  Attenuation  Measurement 

The  B-mode  imager  will  provide  two  important  functions:  a)  acquisition  of  2D  B- 
mode  images  for  construction  of  a  3D  geometric  model  of  patient  anatomy  and  b)  provide 
an  estimate  of  the  ultrasound  absorption  field  in  the  body  tissue  to  provide  starting  point 
parameters  for  the  models. 

Geometric  Model  Construction 

Figure  2  diagrams  the  ultrasound  based  data  acquisition  system  which  consists  of  a 
B-mode  ultrasound  imager  to  acquire  patient  geometry.  This  image  will  serve  as  the 
diagnostic  medical  image  that  the  geometric  model  will  be  based  upon. 

The  construction  of  the  geometric  model  is  the  most  labor  intensive  step  in  solving 
for  the  propagation  of  ultrasound  through  tissue.  In  this  step,  it  is  necessary  to  determine 
the  exterior  boundary,  the  boundary  conditions,  and  the  internal  spatial  distribution  of 
material  properties  and  enter  these  data  into  the  computer  for  analysis.  The  most  common 
technique  in  use  today  is  to  base  the  geometric  model  upon  a  standard  diagnostic  medical 
image  such  as  CT  or  MRI.[5]  Regions  of  interest  are  then  traced  out  over  this  image  and 
the  interior  region  is  meshed  for  the  appropriate  numerical  solution. 
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In  the  first  stage  of  this  work  only  two-dimensional  models  will  be  constructed 
based  upon  a  single  image  from  the  B-mode  scanner  (see  Figure  3).  The  imager  is 
equipped  with  a  software  measurement  package  and  a  hard  copy  output  device.  The 
operator  will  capture  the  image  of  interest,  determine  regions  of  constant  material 
properties,  measure  the  dimensions  of  these  regions,  and  print-out  that  image.  Once  the 
geometric  information  is  obtained,  a  two  dimensional  unstructured  mesh  generator[6]  is 
used  to  mesh  the  domain  with  quadrilateral  elements  for  a  finite  element  analysis.  (Note: 
this  is  all  that  was  planned  in  the  original  proposal.) 


The  next  step  is  to  acquire  a  series  of  spatially  distributed  images  and  assemble 
them  into  a  3D  model.  The  clinical  solution  will  require  the  physician  or  assistant  to  sketch 
the  important  anatomical  features  (skin,  tumor,  bone)  on  a  series  of  2D  “slices”.  The 


input  to  the  model. 


27-11-96  1 

63:36  PM  ■ 

Linear  ■ 

7.5  !1Hz  1 

IdnfllnljEB  ■ 

Figure  3:  This  is  a  typical  image  produced  by  the  B-mode  ultrasound  scanner. 


Ultrasonic  Attenuation  Measurement 

The  power  deposited  at  a  specific  point  in  tissue  by  the  ultrasound  transducer 
depends  on  the  intensity  of  the  sound  field  at  that  point  and  the  absorption  coefficient  of  the 
tissue  at  that  point.  Clearly,  the  intensity  at  a  point  depends  on  the  initial  intensity  and  the 
attenuation  of  the  signal  in  the  tissue.  This  attenuation  is  due  to  scattering,  absorption,  and 
geometric  attenuation.  Previous  research  has  shown  that  the  attenuation  is  highly  variable 
within  human  body  tissue,  with  measurements  varying  significantly  for  different  tissue 
types,  for  similar  tissue  on  different  subjects,  and  even  for  the  same  tissue  on  the  same 
patient  on  different  days. [7, 8, 9,10]  Thus,  it  is  important  to  measure  the  attenuation  of  the 
ultrasonic  signal  in  vivo  on  the  patient  at  the  treatment  site. 

Using  the  data  acquisition  system  shown  in  Fig.  2,  it  is  possible  to  estimate  the 
attenuation  of  the  ultrasound  at  different  points  in  the  patient  tissue  from  the  RF  signal 
collected.  This  is  accomplished  using  a  Log  Spectral  Difference  technique.[8,9,l  1]  Details 
of  this  technique  are  discussed  in  Appendix  A.  It  is  important  to  note  that  Log  Spectral 
Difference  is  one  of  several  possible  techniques  and  provides  an  estimate  of  the 
attenuation.  [6]  For  this  reason  it  is  important  to  test  the  accuracy  of  the  technique  with 
other  measurements. 
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Tests  were  conducted  where  the  acoustic  attenuation  of  a  series  of  samples  were 
measured  using  both  log  spectral  difference  and  the  Force  Balance  technique.  Using  the 
Force  Balance  technique  the  force  generated  when  an  acoustic  wave  impinges  on  an 
absorbing  material  is  measured  and  related  to  the  power  of  the  wave.  Using  this  technique 
it  is  possible  to  make  a  direct  measurement  of  the  attenuation  or  absorbed  power.  For  this 
reason,  the  force  balance  technique  is  used  as  a  standard  to  which  the  log  spectral 
difference  measurements  are  compared.  While  it  is  generally  considered  accurate  and 
reliable,  the  force  balance  technique  is  not  suitable  for  use  in  the  clinical  experiments  due  to 
geometric  constraints. 

Using  the  techniques  described  above,  initial  tests  have  been  run  on  chicken  breast 
to  compare  the  results  of  the  two  techniques  discussed  above.  Using  the  Log  Spectral 
Difference  Technique,  the  normalized  attenuation  coefficient  for  the  chicken  was  measured 
as  0.1608  Nepers  x  cm'1  x  MHz1. [4]  Using  the  Force  Balance  system  measured  values 
ranged  from  0.108  to  0.138  Nepers  x  cm‘x  MHz1.  These  values  are  of  the  same 
magnitude  as  values  reported  in  the  literature  (0.06  -  0.13  Nepers  x  cm'1  x  MHz1). [12] 


Subsystem  2:  Absorbed  Power  Measurement 

While  the  log  spectral  difference  technique  is  convenient  because  it  is  non-invasive 
and  allows  measurement  of  the  attenuation  field  within  the  tissue,  it  does  not  give  exact 
measurements.  Attenuation  field  data  received  from  the  B-scan  system  will  be  used  as  a 
starting  point  for  the  parameters  used  in  the  math  models.  If  the  inverse  optimization 
scheme  is  to  be  successful,  it  is  also  necessary  to  make  precise  experimental  measurements 
of  the  power  deposited  in  the  tissue  by  the  treatment  transducer. 

This  research  is  based  upon  an  existing  clinical  system  which  provides  all  the  tools 
necessary  to  perform  power  deposition  measurements  using  a  thermal  technique.  Details  of 
the  original  technique  can  be  found  in  [13,14]  with  more  recent  applications  discussed  in 
[15,16,17,18].  The  relation  between  the  absorbed  power  due  to  relaxation  (as  measured 
by  the  thermal  techniques)  and  the  acoustic  pressure  calculated  by  the  numerical  models 
was  put  forth  by  Nyborg.[19] 


Figure  4:  This  is  a  schematic  diagram  of  the  operation  of  the  absorbed  power  measurement 
test. 


9 


*  A 


Conceptually  these  tests  will  use  the  schematic  shown  in  Fig.  4.  The  treatment 
transducer  will  be  scanned  at  low  power  to  identify  the  locations  of  the  imbedded 
thermocouples.  When  the  focus  of  the  transducer  is  brought  into  the  vicinity  of  a 
thermocouple,  the  thermocouple  will  exhibit  a  sharp  temperature  rise.  Once  the 
thermocouple  is  located,  the  transducer  focus  will  be  centered  on  that  location,  and  power 
applied.  For  a  given  input  power,  the  rate  of  increase  in  temperature  measured  by  the 
thermocouple  is  a  direct  measure  of  the  absorbed  power  at  that  point  (and  also  related  to  the 
attenuation  in  the  tissue  between  the  transducer  and  the  measurement  point).  The 
temperature  changes  experienced  during  these  tests  should  be  brief,  with  maximum 
temperatures  significantly  below  treatment  levels. 


Subsystem  3:  Ultrasound  Propagation  Modeling 

Results  of  the  two  models  currently  under  investigation  are  discussed  here.  Details 
of  their  development  are  included  in  Appendix  B. 

A  hybrid  numerical  model  using  integral  methods  to  solve  for  the  acoustic 
displacements  in  the  region  where  the  ultrasound  propagates  from  the  transducer  to  the 
patient’s  skin  and  using  a  finite  element  solution  of  the  wave  equation  for  propagation  in 
the  body  has  been  developed  and  preliminary  tests  have  been  performed.  Since  the  model 
predicts  the  displacements  and  stresses  propagating  through  the  tissue,  it  automatically 
accounts  for  reflections  and  transmission  at  interfaces.  As  a  result  we  see  in  preliminary 
results  a  tendency  for  standing  waves  to  be  formed  caused  by  the  artificial  boundary 
conditions  set  at  the  edges  of  the  model  space  (see  Appendix  B  for  details).  The  model 
currently  has  difficulty  handling  the  frequencies  produced  by  the  treatment  transducers. 
This  results  from  the  fact  that  the  finite  element  mesh  must  use  smaller  elements  at  higher 
frequencies.  This  results  in  more  elements  for  the  same  physical  model  size  and  an 
exponential  increase  in  computation  time. 

A  second  model,  based  on  the  parabolic  equation  method  has  also  been  developed 
and  has  undergone  preliminary  tests.  The  parabolic  equation  method  is  a  very  efficient 
method  for  modeling  acoustic  wave  propagation  through  low  contrast  acoustic  materials. 
This  solution  technique  is  somewhat  more  of  an  approximation  than  the  hybrid  model,  but 
because  of  the  source/receiver  geometry  used  in  the  hyperthermia  experiment,  and  the 
relatively  low  contrast  of  breast  tissue,  it  is  possible  to  utilize  this  efficient  approximation. 
The  resulting  speed-up  relative  to  the  finite  element  method  is  dependent  upon  the  contrast 
but  is  found  to  be  fairly  large.  This  is  especially  true  of  the  large  3D  problems  simulated  in 
this  report. 

Speed  Issues: 

We  have  found  that  there  is  difficulty  in  using  standard  mesh-generating  code  to 
model  a  7.5  by  7.5  by  7.5  cm  piece  of  tissue  at  500  KHz.  However,  we  have  carried  out 
simulations  at  the  full  2  MHz  using  the  Parabolic  Fourier  Marching  method.  The  operating 
frequency  of  the  device  is  2  MHz  so  we  are  confident  that  the  present  code  will  serve  as  a 
suitable  model,  provided,  of  course  that  the  accuracy  can  be  verified. 

This  difference  in  speed  is  due  in  part  to  variations  in  the  density  of  spatial  sampled 
required  for  the  various  models.  The  parabolic  marching  method  requires  spatial  sampling 

at  only  X/2,  where  X  is  the  wavelength  in  the  slowest  material  present,  at  the  highest 

frequency.  This  yields  a  savings  of  a  factor  of  23  =  8  over  the  integral  equation  method 
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which  requires  a  sampling  of  X/4.  The  finite  element  solution  requires  a  sampling  of 
approximately  X/8,  so  the  parabolic  method  gives  an  improvement  of  roughly  4^=  64. 

Figures  5  and  6  demonstrate  the  kind  of  information  available  from  the  models. 
They  have  been  duplicated  here  from  Appendix  B  for  convenience.  In  Figure  5  we  see  a 
continuous  plane  wave  propagating  from  the  left,  entering  the  body  through  a  vertical  skin 
layer.  The  wave  travels  through  an  elliptical  fatty  region  and  into  the  spherical  tumor.  The 
left  hand  figure  shows  the  wave  energy  density  that  tends  to  get  weaker  as  the  wave  moves 
to  the  right.  Note  that  the  tumor  absorbs  most  of  the  wave  energy  resulting  in  a  “dark” 
shadow  behind  the  tumor.  The  right  hand  figure  shows  the  energy  absorbed.  Note  that  in 
the  water  and  normal  tissue  the  absorbed  energy  is  low;  even  when  the  wave  energy  is 
high.  Then  there  is  an  increase  in  absorbed  energy  as  the  wave  passes  into  the  tumor.  This 
results  from  the  higher  absorption  coefficient  in  the  tumor.  This  effect  is  shown  more 
clearly  in  Fig.  6  which  is  a  line  plot  of  the  absorbed  energy  along  a  central  axis  in  Fig.  5. 


Figure  5:  Shown  on  the  left  is  the  predicted  wave  energy  density  produced  by  an 
incident  continuous  plane  wave.  On  the  right  is  shown  the  density  of  the 
energy  absorbed  by  the  tissue. 


Model  Accuracy 

While  the  finite  element  model  is  apparently  quite  accurate  it  is  necessary  to 
establish  conditions  on  boundaries  of  the  modeled  region.  Since  these  boundaries  are 
artificially  set  in  the  model,  the  conditions  at  the  boundary  must  be  for  free  space  radiation. 
These  conditions,  as  established  in  our  model,  were  not  sufficiently  accurate,  and  there 
were  standing  waves  set  up  in  the  computational  grid  that  were  spurious.  Changing  model 
geometry  and  boundary  conditions  would  affect,  but  not  eliminate  these  spurious  waves. 
This  problem  with  the  boundary  conditions,  along  with  the  speed  issues  discussed  above 
have  resulted  in  the  rejection  of  the  FE  model  as  an  option  for  our  current  research. 
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Figure  6:  Thermal  energy  deposition  along  a  central  axis  predicted  by  the  parabolic  model. 

The  verification  of  the  parabolic  marching  method  is  discussed  below.  The  "gold 
standard"  in  the  initial  tests  has  been  analytic  spherical  Bessel  function  solutions  for  simple 
geometries  that  have  analytic  solutions.  For  the  complicated  cases  encountered  in  the 
laboratory,  we  will  use  the  integral  equation  solutions. 

The  following  figures  demonstrate  the  agreement  between  the  integral  equation 
method  and  the  parabolic  method.  Figure  7  shows  the  geometry  used  for  this  analysis. 
Figure  8  shows  a  color  mapped  representation  of  a  slice  of  the  field.  Figure  9a  and  9b 
show  the  results  using  a  speed  of  sound  in  the  sphere  significantly  less  than  and 
significantly  greater  than  that  of  water.  Given  the  differences  in  the  two  models,  one 
would  expect  error  to  increase  as  the  mismatch  in  sound  velocity  between  the  tumor  and 
water  increases.  Thus,  these  two  conditions  should  represent  worst  case  errors.  Figure  10 
is  an  axial  plot  of  the  same  data. 


Figure  7 :  Model  geometry  for  accuracy  comparisons. 
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Figure  8:  This  is  a  color  mapped  image  of  the  power  deposition  in  the  model  shown  in 
Figure  7. 


Figure  9:  a)  The  speed  of  sound  internal  to  the  sphere  was  1267.7  m/sec,  substantially 
less  than  that  which  would  be  encountered  in  the  breast,  thereby  making  the 
impedance  contrast  quite  a  bit  larger  than  is  likely  to  be  encountered  in  practice, 
b)  Here  the  speed  of  sound  internal  to  the  sphere  was  1677  m/sec,  which  is 
substantially  greater  than  any  tumor  that  would  be  found  in  the  breast.  Here  the 
solid  line  represents  the  integral  model  and  the  dashed  line  represents  the 
parabolic  model. 


Comparison  of  axial  field 


Figure  10:  The  integral  equation  is  shown  in  solid,  and  shows  the  diffraction  pattern  due 
to  reflections  from  the  hard  sphere.  This  plot  shows  quantitatively  that  the 
parabolic  method  does  not  take  into  account  the  reflections  of  sound  waves. 
Conditions  for  this  graph  correspond  to  those  in  Fig.  9b. 

The  excellent  match  between  the  two  models  in  the  forward  scattered  diffraction 
patterns  shown  in  Fig.  9  and  the  poor  match  in  the  reflected  diffraction  pattern  shown  in 
Fig.  10  highlights  the  differences  between  the  two  models.  The  excellent  overall  match 
between  the  two  models  in  this  worst  case  condition,  coupled  with  the  enormous  speed 
enhancement  associated  with  the  parabolic  model  makes  it  an  excellent  candidate  for  our 
application. 


Subsystem  4:  Inverse  Optimization  Technique 

Inverse  Techniques 

Inverse  problems  arise  when  the  physical  environment  either  severely  limits  or 
completely  prohibits  the  direct  measurement  of  the  parameters  necessary  to  uniquely 
determine  a  system's  dependent  variables.  When  using  ultrasound  hyperthermia  to  treat 
breast  cancer 

the  parameters  of  interest  are  the  spatial  distribution  and  values  of:  (a)  tissue  density,  (b) 
acoustic  pressure  amplitude  attenuation  coefficient  and  (c)  wave  speed.  The  wave  speed 
and  acoustic  pressure  amplitude  attenuation  coefficient  can  be  treated  as  one  quantity 
in  terms  of  the  complex  wave  number.  The  dependent  variable  of  interest  is  the  absorbed 
power,  (q). 

If  one  knows  all  the  parameters  of  a  problem  then  it  is  possible  to  solve  the 
"forward"  problem  for  the  dependent  variables  using  the  differential  wave  equation.  The 
inverse  problem  arises  when  one  does  not  know  the  parameters  but  is  able  to  measure  some 
other  property  of  the  system,  usually  some  of  the  dependent  variables.  For  ultrasound 
hyperthermia  of  breast  cancer  it  is  possible  to  pre-determine  geometric  regions  of  constant 
property  (see  Subsystem  1)  and  then  measure  on-line  the  value  of  the  absorbed  power  (see 
Subsystem  2)  at  a  limited  number  of  internal  points  —  the  thermocouple  junctions. 
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The  inverse  technique  developed  for  this  research  utilizes  a  standard  inverse  technique  for 
non-linear  problems[20]  that  is  based  on  minimizing  the  error  between  the  model  estimated 
absorbed  powers  and  the  clinically  measured  absorbed  powers.  The  minimization  is 
accomplished  by  changing  the  value  of  the  unknown  parameters  in  the 
model  until  the  “cost"  function, 


Nm 

O(e)  =  ^  ef  Equation  [1] 

i=l 

where 

ei={Qv)i  “Wv/j 

is  minimized.  This  technique  is  also  termed  parameter  estimation.  The  inverse  process 
begins  by  assuming  some  value  for  the  unknown  parameters  and  solving  the  “forward" 
problem.  The  error  between  the  model  estimated  dependent  variables  and  the  measured 
dependent  variables  and  the  resulting  changes  in  absorbed  power  are 
given  by  the  following  linear  approximation: 

[J]  (Ay)  =  (e)  Equation  [2] 

The  “Jacobian"  matrix,  J,  is  constructed  from  numerical  simulations  and  is  given  by 


dMi 

Sy: 

dy2 

dyNp 

4?v)2 

dM2 

dy2 

dyNp 

Nm 

dMNm 

dMNm 

dyl 

dy2 

- 1 

£ 

Equation  [3] 


The  inverse  method  utilized  herein  follows  the  following  cycle  of  calculations  in  which  the 

model  parameters,  y,  are  updated  according  to  equation  [2]  until  there  are  negligible 
changes. 

1.  Using  the  current  values  of  model  parameters,  y,  calculate  the  Jacobian  in  equation  [3] 

2.  Construct  the  error  vector  using  the  current  values  of  model  parameters,  y 

3.  Invert  equation  [2]  using  singular  value  decomposition  (SVD)  to  solve  for  changes  in 
model  parameters,  y 

4.  Update  model  parameters 

5.  Stop  if  the  change  in  model  parameters  is  small  or  the  error,  e,  is  acceptable 
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Future  Work 

Third  year  efforts  will  involve  parallel  efforts  in  all  four  subsystems: 

Subsystem  1:  A  hardware  system  will  be  developed  to  allow  the  geometry  acquisition  to  be 
automated.  Software  development  will  continue,  moving  toward  defined 
clinical  methods  for  geometry  acquisition  and  attenuation  estimation. 

Subsystem  2:  Software  for  the  experimental  plan  discussed  above  will  be  implemented. 
Clinical  tests  for  this  subsystem  are  anticipated. 

Subsystem  3:  A  final  model  will  be  selected  and  work  will  begin  to  produce  an  automatic 
interface  between  Subsystem  3  (the  model)  and  the  other  subsystems. 

Subsystem  4:  Software  development  will  continue.  Beginning  with  simulated  data,  then 
using  clinical  data  when  it  is  available  from  Subsystem  2,  the  algorithm  will 
be  debugged  and  tested. 
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Appendix  A:  LOG  SPECTRAL  DIFFERENCE 
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The  log  spectral  difference  method  is  a  technique  of  determining  an  ultrasound 
pressure  amplitude  attenuation  coefficient  (UPAAC),  based  on  the  power  spectrum  of  a 
reflected  ultrasound  signal.  [8,9,1 1]  The  idea  behind  this  technique  is  to  compute  the 

normalized  attenuation  coefficient,  [3,  by  determining  the  reduction  in  amplitude  of  the 

spectrum.  The  UPAAC,  a,  is  then  calculated  from  Equation  Al,  where  (3  is  the  attenuation 
coefficient  normalized  to  both  frequency  and  depth,  and  f  is  the  frequency.  Depending  on 
the  type  of  tissue,  n  has  been  found  to  be  between  1  and  1.3. [21] 

oc  =  /3f n  Equation  [Al] 

To  simplify  the  calculations,  a  linear  relationship  between  a  and  f  (n=l)  is  assumed.  The 
process  of  calculating  (3  will  be  discussed  in  the  following  paragraphs. 

Determining  (3  begins  with  analyzing  the  reflected  radio  frequency  (RF)  signal  of  a 
pulse  sent  from  an  ultrasound  transducer  into  the  tissue  of  concern.  Portions  of  the  signal 
originating  from  the  front  and  back  areas  of  the  region  of  interest  are  chosen.  The  shallow 
region  is  selected  to  be  the  tissue  segment  just  after  a  tissue  interface,  because  the  irregular 
geometry  and  surface  roughness  at  the  interface  will  scatter  the  ultrasound  wave  and 
produce  inaccurate  results.  The  deep  region  is  chosen  before  a  tissue  interface  for  the  same 
reason.  The  optimal  distance  between  the  two  chosen  regions  has  been  shown  by  Kuc  to 
be  two- thirds  of  the  region  of  interest  (see  Figure  A 1).  [22,23]  A  typical  RF  signal  received 
from  a  pulse  sent  into  tissue  is  shown  in  Figure  A2.  The  portions  of  the  signal  associated 
with  the  shallow  and  deep  regions  are  marked  on  the  figure.  However,  these  segments  of 
the  signal  are  not  isolated. 


Figure  Al :  The  ellipse  represents  a  tissue  region  of  interest  (tumor,  organ,  etc.). 
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Backscatter  Radio-Frequency 


Data  Points 


Figure  A2:  This  represents  the  RF  waveform  of  a  single  A-line  collected  from  the  B-scan 
imager.  Two  regions  selected  from  the  object  in  Fig.  A1  are  indicated. 


The  regions  selected  from  the  RF  signal  contain  overlapping  reflections  from 
different  areas  of  tissue.  To  reduce  the  many  complexities  associated  with  this 
phenomenon,  the  signal  from  each  region  is  divided  into  overlapping  sections,  and  each 
section  is  windowed  using  a  Hanning  or  similar  window.  Next,  the  log  of  the  Discrete 
Fourier  Transform  (DFT)  is  calculated  for  each  section  to  form  the  power  spectrum.  This 
spectrum,  exhibited  in  Figure  A3,  is  for  the  shallow  and  deep  segments  of  Figure  A2. 

In  the  power  spectrum  form,  the  signal  can  be  quickly  analyzed.  The  power 
spectrum  of  the  deep  region  is  subtracted  from  the  spectrum  of  the  shallow  region,  leaving 
a  single  line  representing  the  log  spectral  difference  (thus  the  name  of  the  method).  The 

slope  of  this  line,  co,  is  calculated  using  linear  regression.  Figure  A4  shows  the  result  after 
the  power  spectrum  of  the  deep  and  shallow  segments  have  been  subtracted  and  linear 
regression  has  been  performed.  After  the  slope  of  the  line  is  found,  it  is  divided  by  the 
round-trip  distance  between  the  chosen  segments  giving  the  normalized  attenuation 

coefficient,  P,  as  shown  in  Equation  A2. 


Equation  [A2] 


The  variables  X|  and  x2  represent  the  distances  to  the  shallow  and  deep  segments 

respectively.  After  P  is  determined,  the  UPAAC  can  be  calculated  as  demonstrated 
previously. 
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Figure  A3:  The  log  of  the  magnitude  of  the  RF  waveform  for  the  deep  and  shallow 
regions  is  plotted  as  a  function  of  frequency. 


Linear  Fit  to  Loss 


Figure  A4:  By  subtracting  the  two  curves  from  Figure  A3,  the  Log  Spectral  Difference  is 
formed. 

The  UPAAC  calculated  from  a  single  RF  signal  is  not  a  reliable  measurement  of 
acoustic  attenuation.  Discontinuities  and  other  distortions  in  the  tissue  will  skew  the  results 

if  only  one  RF  signal  is  analyzed  to  determine  a.  To  avoid  any  anomalies  that  may  be 
present  in  a  single  signal,  the  results  from  the  analyses  of  several  contiguous  signals  should 
be  averaged  together.  This  will  produce  a  more  accurate  and  precise  method  of  determining 
the  UPAAC  value  for  a  given  substance. 
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Appendix  B:  ULTRASOUND  PROPAGATION  MODELS 

Parabolic  Equation  Fourier  Marching  Method 

Introduction 

The  parabolic  equation  method  is  a  very  efficient  method  for  modeling  acoustic 
wave  propagation  through  low  contrast  acoustic  materials  (such  as  breast  tissue).  The 
original  or  classical  method  requires  for  its  applicability  that  energy  propagate  within 
approximately  ±20°  from  the  incident  field  direction.  Later  versions  allowed  propagation  at 
angles  up  to  ±90°  from  the  incident  field  direction.  Further  modifications  provide  accurate 
backscattering  information,  and  thus  are  applicable  to  the  higher  contrasts  encountered  in 
nondestructive  imaging,  EM  and  seismic  applications.  [24,25] 

The  source/receiver  geometry  used  in  the  hyperthermia  experiment,  and  the 
relatively  low  contrast  of  breast  tissue,  allow  us  to  utilize  this  efficient  approximation.  The 
resulting  speed-up  relative  to  the  finite  element  method  is  dependent  upon  the  contrast  but  is 
found  to  be  fairly  large.  This  is  especially  true  of  the  large  3D  problems  simulated  in  this 
report.  We  have  found  that  there  is  difficulty  in  using  standard  mesh-generating  code  to 
model  a  7.5  by  7.5  by  7.5  cm  piece  of  tissue  at  500  KHz.  However,  we  have  carried  out 
simulations  at  the  full  2  MHz  using  the  Parabolic  Fourier  Marching  method.  The  operating 
frequency  of  the  device  is  2  MHz  so  we  are  confident  that  the  present  code  will  serve  as  a 
suitable  model,  provided,  of  course  that  the  accuracy  can  be  verified.  The  verification  of 
the  parabolic  marching  method  has  been  successfully  carried  out  for  smaller  problems,  but 
is  still  be  tested  for  the  larger  3D  problems  that  are  required  for  this  report.  The  "gold 
standard"  in  the  initial  tests  has  been  analytic  spherical  Bessel  function  solutions  for  simple 
geometries  that  have  analytic  solutions.  For  the  complicated  cases  encountered  in  the 
laboratory,  we  will  use  the  integral  equation  solutions.  Furthermore  the  coarse  grain 
parallelization  employed  in  the  integral  equation  method  is  equally  applicable  to  the 
parabolic  algorithm. 

The  Parabolic  equation  approximation,  and  more  exactly,  the  "split  step  Fourier 
Method"  is  a  new  method  based  upon  earlier  literature.  [26,27] 

The  source/receiver  geometry  requirements  and  the  relatively  low  contrast  of  breast 
tissue,  allow  us  to  utilize  this  efficient  approximation,  for  breast  cancer  scanner  and 
hypothermia  applications,  for  example.  In  particular  because  we  may  elect  to  use  only 
transmission  data  in  the  breast  hyperthermia  problem,  we  are  able  to  use  this  "parabolic 
approximation"  in  a  straight-forward,  simple  manner  [26]. 

The  Parabolic  method  can  be  adjusted,  however,  to  include  the  effects  of  backscatter, 
although  this  approximation  increases  the  computational  complexity  of  the  algorithm. 


Discussion  and  Derivation  of  Fourier  marching  method 

To  derive  and  elucidate  the  parabolic  equation  method,  we  begin  with  the  2-D 
Helmholtz  wave  equation  governing  wave  propagation  in  inhomogeneous  media: 


( 


dx‘ 


+ 


k2{x,y)  + 


\2  \ 


\f(x,y)  =  0 


Equation  [Bl] 
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Now  Fourier  transforming  y  to  X  results  in: 


1  A 


-y+  -—kza)*-K  \\f(X,X)  =  0 


where: 


dxJ  V2n 


k2  (\)  =  k2  (x,X)  =  \k2(x,y)e  lXydy 


Equation  [B2] 


Equation  [B3] 


and  the  *  notation  denotes  convolution: 


f(X)*g(X)=  if(x-x')g(X')dk' 


Eqn.  (B2)  can  be  factored  in  the  sense  of  pseudo-differential  operators  [28] 


Equation  [B4] 


d  I!  'S  7Y  d  n  7  - 

- +  L — k2(X)*-X2 - L — k2  (X)*-X2  f(x,X)  =  0 

dx  V  2n  J  a*  V  2jc  f 


Equation  [B5] 


An  intuitive  feel  for  the  manner  in  which  the  Parabolic  equation  arises  can  be  seen  by 
looking  at  particularly  simple  case:  when  k  is  a  constant,  then  we  have,  upon  Fourier 
Transforming  (Bl): 


iv k2  -  A2  \f(x,A)  =  0  Equation  [B6] 


where  the  square  root  is  now  simply  the  square  root  of  a  scalar.  From  B6,  it  is  clear  that  in 
this  special  case,  two  solutions  exist: 


-^—  +  Wk2-; l2  f(x,X)~  0,  f(x,X)  =  g(A)e  r  Equation  [B7] 

\dx  J 

for  an  arbitrary  function  g,  which  represents  a  right  moving  wave  for  eim  time  dependence 
and: 

(4 —  i^Jk2-X2\f(x,  X )  =  0,  f{x,  A )  =  g(A  Equation  [B8] 


- - Hk  -Af  \fix,X)  =  0,  f{x,X)  =  g{X)el> 


representing  a  wave  which  "moves"  to  the  left.  Suppose  that  we  know  that  the  field  is  due 
to  sources  entirely  on  the  LHS  of  the  x-y  plane.  Then,  for  x>0  we  must  have: 


f{x,X)  =  g{X)e~ 


Equation  [B9] 


Note  that  knowledge  of  f  on  the  line  x=0  (boundary  condition)  completes  the  solution 
since: 

/( 0,  A)  =  g(X)  Equation  [B 1 0] 
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and  so,  on  inverse  Fourier  Transforming: 

■4  00  _ 

f(x, y)  =  —  f /(0,  X  )e~ix4^eiXydX ,  x  >  0  Equation  [B 1 1] 
27T 

In  fact,  for  any  xo>0: 

a  _ 

/( x,y)  =  —  J /(jc0, X)e~iix~Xa)4iF^ eiXydX,  x>x0>0  Equation  [B12] 

2  JT 

— oa 

The  idea  behind  the  parabolic  equation  method  is  to  try  to  factor  a  general  inhomogeneous 
(i.e.  k=k(x,y))  problem  into  forward  and  backward  moving  (in  x)  factors  and  then  solve 
only  the  forward  (+x)  moving  part  (assuming  that  the  field  source  is  to  the  left  (on  the  x- 
axis)  of  the  scatterer). 

Let  xn  =  nA ,  n=0,...  then  B 12  is: 

fn+l(} 0  =  ~  ]  fna)e~i&^eiXydX  Equation  [B13] 

Since  B13  is  local  to  x=(n+l/2 )A  we  consider  the  discretization/approximation: 

00  .  I  2  « ? 

fn+i(y)  =  T-  e^dK 

2ft  -°o  Equation  [B14] 

for  propagating  forward  through  k  that  is  inhomogeneous  in  x  and  y  i.e.  k~  k(x,y),  and 
fn(y)~  f  (xn>y )•  Computationally  (B 14)  would  be  an  inverse  Fourier  transform 
except  for  the  y-dependence  of  k. 

Defining: 

k2n  =k2  (xn+ ll2) 

gives: 

/  00  •  172  v2 

f»*l(y)  =  2-  \fnO-)e  ^  eA,dX 


Advanced  Parabolic  Marching  Method 

Now  we  desire  to  take  the  y  -  dependence  out  from  under  the  square  root  in  order 
that  it  may  then  be  factored  out  from  under  the  integral,  which  will  result  in  the  integral 
being  an  inverse  Fourier  Transform,  and  yield  a  substantial  increase  in  efficiency.  This 
goal  can  be  achieved  approximately  in  the  following  manner: 
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e  lAkn,m^lJ  lkn,m )  _  g  ^(fn.m  kn  ) ^ ~iAkn -jl-{h/kn f 


so  that,  defining  kn  m  =  k(xn+1/2  ,  ym  )  *  and  Vm  =mA>  then  yields: 

;  J/„  (X ea*»dX 


fn+1  (ym) 


2n 


or 


-iA  £ 


cn,m  kn) 


fn+1  (  y  m  ) 


2n 


J/„  (X  )e~iMnf^y2  e*ym  dX 


Equation  [B15] 


We  have  thus  approximately  maintained  the  y  variation  in  k,  while  simultaneously 
giving  a  Fourier  Transform  formulation  for  the  field  fn+j. 


riA(kn,m-kn>  «, 


fn+i(y) 


iMc^l-\ \ 


2n 


I  fn(X)e  "V  t/+J  ea.,dX 


and  in  fact  using  the  notation: 


f(y)  =  ^~  jfa)e-iXydi 

2tc_00 


to  indicate  the  Inverse  Fourier  Transform,  gives: 

A 

fn+l(y)  =  fn(^)Pn 

where  Pn  is  the  range  dependent  propagator  defined  as: 


P  (-y)  =  e-i^k2(xn+u2)-'k2 

This  is  one  basic  equation  for  the  parabolic  split  step  method.  Various  alternative  forms 
can  be  used  in  similar  algorithms.  A  common  form  for  the  parabolic  equation  is  derived  by 
using  the  binomial  approximation: 


in  the  above  integral.  This  yields  a  "standard"  form  of  the  split-step  parabolic  equation 
method: 
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a)e  2heay”dX 

Equation  [B16] 

Numerical  experiments  have  indicated  the  superiority  of  B15  over  B 16.  In  Fourier 
notation,  B16  is: 


fn+I(ym)  =  e~‘Un'mF-1\ 


iA- 


e  2k°F{fn}CK)\(ym) 


Equation  [B17] 


This  is  the  usual  form  of  the  split-step  PE  method.  Our  more  general  equation  is  (in 
Fourier  notation): 


/„+,  (ym)  =  F~'  {  )(X) 


<ym) 


An  interesting  interpretation  of  the  split  step  method  is  evident  from  this:  The  algorithm 
step  consists  of  an  exact  propagation  a  distance  of  A  in  kQ  (which  includes  diffraction) 
followed  by  a  phase  shift  correction,  i.e.,  multiplication  by 

e~iA{k^~h)  Equation  [B 18] 

which  corrects  the  phase  of  a  plane  wave  traveling  forward. 


ISUtniferical  Experiments;. 


The  geometry  for  the  first  experiment  consisted  of  Figure  B1  —  which  has 
cylindrical  symmetry  about  the  horizontal  axis.  The  tumor  was  assumed  to  be  perfectly 
spherical  with  a  diameter  of  4  cm.  it  was  placed  on  the  center  axis  with  the  center  of  the 
sphere  located  at  a  distance  of  4  cm  from  the  tissue/water  interface.  The  transducer  was 
assumed  to  be  located  on  the  central  axis,  with  its  center  located  10  cm  from  the 
tissue/water  interface.  The  incident  field  in  the  first  experiment  was  assumed  to  be  a  plane 
wave  propagating  along  the  y  axis.  The  z  axis  is  the  vertical  axis,  the  x-axis  is  coming  out 
of  the  page,  and  the  y-axis  is  the  horizontal  axis. 


Once  the  pressure  field  has  been  calculated  via  the  Fourier  Domain  Parabolic 
Propagation  method,  the  equivalent  internal  thermal  energy  generation  due  to  ultrasound 
absorption  can  be  calculated  using  the  formula: 


Q  =  a 


Z 


The  formula  I  use  for  the  complex  object  function  y  is: 


YfX)=f*!wl  _i=fehi^)T_i 

K  J  l  kCo  ) 


where  the  complex  wavenumber  kc(x)  is  the  spatially  dependent  wavenumber  which 
includes  the  attenuation  coefficient  as  its  imaginary  part.  The  real  part  of  kc(x)  is  the 
standard  wavenumber  co/c(x),  where  c(x)  is  the  spatially  dependent  speed  of  sound. 
Therefore  the  object  function  ytakes  the  form: 


=  |g?/£(x)-t«(x)| 

l  J 


Y(x) = f: ZllWzlZ wf  _1=f_£e__  ,■„(*)£*  f  - 1 
l  ®lc0  J  U(x)  0)1 


0=1476  m/sec 
p=0.9  g/cm 

a=  0.6  nprs/m/Mms 


c=1650  m/sec 
p=1 .1  g/cm 
— 5  nprs/m/MHz 


c=1500  m/sec 

p=1 .0  g/cm 

a=  0.01  nprs/m/MHz 


c=1570  m/sec 
p=1 .01  g/cm 

a=  1 .0  nprs/m 


Figure  B 1 :  Schematic  of  the  geometry  used  in  the  model  calculation. 
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The  following  table  lists  the  values  of  the  attenuation  coefficient  in  Nepers/m,  and  real 
speed  of  sound  in  m/sec. 


Tissue  type 

a  (Nprs/m) 

c  (m/sec) 

fat 

0.6 

1476 

parenchyma 

1.0 

1570 

tumor 

5 

1650 

water 

0.01 

1500 

Plane  Wave  Incident  Field  Simulation; 

The  following  pictures  show  the  results  from  the  simulation  of  a  plane  wave 
impinging  upon  the  tissue/water  interface  (the  first  vertical  line  encountered  moving  from 
left  to  right),  the  fatty  ellipse,  (shown  in  outline  in  this  cross-section),  and  the  spherical 
tumor  (shown  in  cross  section  as  a  circle).  The  color  map  is  shown  below  the  cross- 
section,  the  z  axis  is  vertical,  the  x-axis  comes  out  of  the  page,  and  the  y-axis  goes  left  to 
right  in  the  plane  of  the  paper. 


Figure  B2:  Shown  on  the  left  is  the  predicted  wave  energy  density  produced  by  an 
incident  continuous  plane  wave.  On  the  right  is  shown  the  density  of  the 
energy  absorbed  by  the  tissue. 
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Figure  B3:  Acoustic  field  strength  along  a  central  axis. 


distance  along  axis  (cm) 
Figure  B4:  Thermal  energy  deposition  along  a  central  axis. 


Focused  Incident  Field 

Now  we  look  at  thermal  energy  deposition  in  the  case  of  an  incident  field  which  is  focused 
at  the  center  of  the  spherical  tumor. 
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Figure  B8:  Acoustic  field  strength  along  a  central  axis. 
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Hybrid  Numerical  Method 

The  hybrid  numerical  method  is  a  combination  approach  to  solving  for  the 
pressures  arising  from  the  wave  equation  for  the  propagation  of  sound  waves  through 
heterogeneous  materials.  The  method  uses  the  Kirchoff  Diffraction  Integral  to  calculate  the 
pressures  in  the  homogeneous  water  bath  region  and  the  Finite  Element  Method  (FEM)  for 
the  inhomogeneous  tissue  region  (see  Figure  B9). 


FEM 


Patient 


Water  Bath 

Kirchoff 

Integral 


Ultrasound 

Propagation 


Ultrasound  Transducer 


Figure  B9:  Diagram  of  the  clinical  situation  for  a  hyperthermia  treatment  of  the  breast  or 
chest  wall  area. 

The  ultrasound  wave  equation  for  the  complex  pressures  is  given  by  the  following 
equation: 


V- 


P  (*) 


Vp(x) 


k  (x) 

+  -r-L^p(x)  =  0, 


P*(x)' 


Equation  [B19] 


the  pressure,  p  and  wave  number,  fcare  complex  quantities.  To  solve  this  equation,  a 
hybrid  method  is  used  which  takes  advantage  of  different  numerical  methods  that  are  more 
appropriate  for  different  parts  of  the  domain. 


Kirchoff  Diffraction  Integral: 

For  the  homogeneous  water  bath,  the  Kirchoff  Diffraction  Integral  (B20)  is  used  to 


1  JkR 

P(x)  =  -T-|-TTn/- 


4 n{  R 


Vp  +  ik 


{  i 


1  +  ^ 

V  kRj 


ip 


Ida' 


Equation  [B20] 
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solve  for  the  pressures  arising  from  the  ultrasound  transducer.  The  primary  purpose  of  this 
integral  solution  is  to  generate  boundary  conditions  at  the  water  bath  -  membrane  interface 
that  are  used  by  the  Finite  Element  Method  described  below.  The  following  assumptions 
are  made  in  this  formulation:  no  reflections  and  no  refractions  at  the  water  bath  -  membrane 
interface.  This  method  is  fast  and  changes  in  the  tissue  regions  do  not  impact  the  solution 
of  (B20).  Consequently,  only  a  smaller  region  of  the  overall  domain  needs  to  be 
considered  when  material  properties  and  tumor  locations  vary. 

Finite  Element  Method  (FEM) 

The  solution  of  the  acoustic  wave  equation  (B19)  can  be  achieved  by  directly 
solving  the  governing  partial  differential  equation  through  application  of  appropriate 
numerical  techniques  such  as  the  finite  element  method.  There  are  typically  two  stages  of 
this  method,  the  preprocessing  stage  and  the  solution  phase.  The  preprocessing  stage  is 
accomplished  using  an  unstructured  mesh  generator,  FASTQ,  to  generate  a  set  of  non¬ 
overlapping  quadrilateral  elements.  The  solution  phase  can  be  developed  utilizing  any 
standard  finite  element  method  since  the  wave  equation  is  self-adjoint  [29,30].  In  this  case, 
application  of  the  methods  of  weighted  residuals  results  in  the  following  identity. 


f  1 

VP*(X) 


Vp(x)J 


+ 


P(x) 

p*(x) 


p(x) 


m=o. 


Equation  [B21] 


Integrating  by  parts  produces  the  “weak  statement”  of  the  previous  identity, 


| V  VVcf ~J~  Vp(x)l  +  p(x)  dQ=  [W—VpdQ  Equation  [B22] 

J  p  (x)  n  'v'  J  n 


p(x) 


Discretizing  p  in  the  usual  finite  element  manner  produces  the  following  linear  algebraic 
relation 


1 


V[W]7 


-j- V[/V]]+[H']r  dn{p(x)}  =  M{p(x)}  =  0. 

P  y  p  [Xj 


Equation  [B23] 


The  complex  wave  equation  does  possess  one  additional  complexity  not  mentioned  in  the 
usual  finite  element  texts  -  the  variables  are  complex  numbers. 
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